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An algorithm and a computer program in Fortran 95 are presented which enumerate the 
Hugenholtz diagram representation of the many-body perturbation series for the ground 
state energy with a two-body interaction. The output is in a form suitable for post- 
processing such as automatic code generation. The results of a particular application, 
generation of IATgX code to draw the diagrams, is shown. 
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1. Many-Body Perturbation Theory 

Many problems of interest in quantum mechanics have no known analytic solution, 
and some approximation is necessary to obtain a solution. A particularly useful 
approximation method is perturbation theory, which consists of finding a way of 
separating the Hamiltonian into a part which a solvable, either numerically or an- 
alytically, and a part which is small (See e.g. 1 for a general treatment). One may 
then start from the solution of the solvable part and build in successive orders of 
the perturbation in a power series approach to obtain a closer approximation to the 
solution of the full problem, providing the series converges. 

In many-body quantum mechanics, the perturbation series is convenently resp- 
resented in diagrammatic form, based on the work of Feynman 2 , Goldstone 3 and 
Hugenholtz 4 . Diagrams may be prescriptively expressed in algebraic form which 
give a complete representation of the perturbation series. One may construct such 
series of diagrams and expressions for different observables in Quantum Mechanics. 
This paper is concerned with the ground state expectation value of the Hamiltonian 
operator that is, the ground state energy, or vacuum amplitude for a many-fermion 
system in the special, but common case that the unperturbed problem is given by 
the solution of the Hartree-Fock equations and the system is governed by a two- 
body interaction. This situation is widely dealt with in textbooks 5:6 . Each order of 
correction in perturbation theory for the ground state energy involves, in general, 
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many diagrams and it is important to be able to enumerate them, order by order, 
if one is to perform a calculation 7 . 

2. Diagrammatic Representation 

The perturbation series for the ground state energy can be given by diagrams in 
the Hugenholtz vein constructed according to the following rules (see 5 for a more 
detailed explanation): 

For an n th order diagram, n dots (vertices) are drawn in a column. These 
vertices are connected by directed lines subject to the conditions 

1. Each vertex has two lines pointing in and two pointing out 

2. Each diagram is connected, i.e. one must be able to go from any one vertex 
to any other by following some number of lines 

3. No line connects a vertex with itself 

4. Each diagram is topologically distinct 

For each diagram, a simple prescription exists to go from the pictorial representation 
to an algebraic expression for the contribution to the total energy from the particular 
diagram 5 . 

Following these rules, one can determine that there are no first order diagrams 
since a line may not link a vertex with itself. There is one second order diagram, 
which is shown in figure 1 and three third order diagrams, which are shown in figure 
2. 




Figure 1: Second order diagram 




Figure 2: Third order diagrams 
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3. Diagram finding algorithm 

A diagram of a given order n or d can be fully described by taking all possible 
(unordered) pairs of distinct vertices, of which there will be 

-zn or d{n ord - 1) (1) 

and specifying the number of lines linking each pair, along with the direction of 
each line. It is clear that the since the number of lines connected to a given vertex 
must be four, that no more than four lines may link any pair. In fact, the case of 
four lines can exist only in the second order diagram shown in figure 1 since any 
two vertices connected by 4 lines can not be connected to any other vertices in the 
diagram, which is inadmissible by point 2 in the above list of rules. Given these 
facts, the following algorithm can be used to find all possible diagrams: 

• Create an ordered list of n pairs numbers with a specification for which vertices 
are associated with each pair: 1&2, 1&3, . . ., l&n or£ 2, 2&3 , 2&4 , . . ., 2&n or( 2, 
• ■ •, (n or d - l)&n 

ora- 

• Allow each number in the list to take on independently the values 0,1,2,3 
corresponding to the number of lines linking the pairs. There will be n^ airs 
such combinations. The case of n or d er — 2 is treated specially by allowing the 
number of lines linking pairs to take the value 4. 

• For each combination of number of lines linking pairs, rule 2 from the above 
list is checked along with the condition that four lines emanate from each 
vertex (due to rule 1). Any combination not satisfying the rules is rejected. 
Any combination which does satisfy these rules is a valid unlabelled diagram. 
The diagram is labelled in all possible ways by adding arrows to the lines in 
all possible permutations of upward- and downward-going arrows and check- 
ing these permutations against rule 1. Those which pass are valid (labelled) 
Hugenholtz diagrams. The diagram is completely specified by the ordered list 
of numbers of lines connecting the pairs, and by a second list which gives the 
number of such lines which are pointing up; the rest must point down. 

This provides an exhaustive search and is guaranteed to find all diagrams, at the 
cost of potential slowness as the number of permutations to check grows rapidly with 
order. The two lists of numbers given in the output may be post-processed to give, 
for example, Fortran code which implements the algebraic form of the diagrams. 
Note that condition 3 is automatically satisfied by the fact that in the specification 
of the pairs, we don't consider pairing a vertex with itself. Similary condition 4 is 
automatically satisfied by the specification of the problem in terms of vertices being 
ordered in a column in space. 

4. Description of Code 
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A Fortran 95 code to calculate valid diagrams using the approach described in 
the previous section has been written and is presented in Appendix A. The code is 
described in terms of the subprograms as follows 

4.1. program hugenholz 

n or d is read in, from which n pa irs is calculated. The pairs structure is set up 
so that the first pair links vertices one and two and so on as described in section 
3. A structure of n pairs numbers is initialized with all numbers zero and passed to 
subroutine new to enumerate the diagrams. 

4.2. subroutine new 

Given a diagram specification, a loop is made over the first number of pairs over 
possible allowed values (0 ... 3 or 4) . The remainder of the diagram specification is 
then passed recursively to the subroutine new for the further numbers to be looped 
over. If we are at the last element of the specification, the diagram is checked for 
consistency using the function consistent and then printed out prepended by '+'. 
Finally a call is made to label to find all ways of labelling the lines with arrows. 

4.3. function consistent 

This function, which returns a logical, takes a possible configuration of lines 
connecting pairs of vertices, and checks first that it satisfies the condition from rule 
1 that each vertex has four lines connected to it. Then the diagram is checked for 
being fully connected by starting at the bottom of the diagram, at vertex 1 and 
following lines until all vertices have been reached, in which case the diagram is 
consistent, or all lines have been followed without reaching some vertices, in which 
case it is inconsistent with rule 2 and we reject it. 

4.4. subroutine label 

A new array is set up which stores the number of upward-pointing lines for each 
pair of vertices. This is initialized to zero and passed to newlabel. 

4.5. subroutine newlabel 

Like new, this subroutine loops over all possible numbers of upward-pointing 
lines that each pair may admit given the number of total lines. The looping is again 
recursive, and once the last pair of lines is reached, the labelled diagram is tested 
with the function testlab. If it passes the number of upward-pointing arrows is 
printed prepended by 

4.6. function testlab 

This function tests for rule 2. The number of lines entering and leaving each 
vertex is calculated. If, for every vertex, both these numbers are two, then the 
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diagram is valid. 

5. Sample output 

The output of the program for second order diagrams is 

+ 4 

* 2 

which says that there is one unlabellcd diagram, which has one pair linked by four 
lines. There is furthermore only one way of labelling it, with two of the lines pointing 
up, as shown in figure 1. The output for n or d = 3 is 
+ 222 

* 020 

* 111 

* 202 

which says that there is again only one unlabellcd diagram, which has three pairs, 
each linked by two lines. This time there are three possible ways of labelling the 
lines with arrows; 020 in which all lines point down except those linking vertices 
1 and 3, 111 in which each pair of vertices is connected by one upgoing and one 
downgoing line, and 202 in which all lines point up except those linking vertices 1 
and 3. These three possibilities are shown in figure 2 from left to right in the order 
given here. 

The graphical representation of the diagrams in this paper were automatically 
generated from the output of the Fortran program hugenholtz by means of a perl 
program* which outputs IATj^X source using the FeynMF package 8 . While it would 
have been quite easy to work out second and third order diagrams by hand, the 
ability to automatically enumerate and represent all diagrams for higher orders 
becomes very useful. For example, the number of (labelled) diagrams at fourth, fifth, 
sixth and seventh orders are 39, 840, 27,300 and 1,232,280 respectively, which would 
be tedious to say the least to enumerate by hand. As an example of a less trivial 
result, the fourth order diagrams, as generated from the output of hugenholtz are 
shown in figure 3 

6. Summary 

An algorithm has been presented and implemented which enumerates and repre- 
sents Hugenholtz diagrams for the ground state perturbation series for a many-body 
system interacting under the influence of a two-body interaction. The representation 
is in a form suitable for post-processing. An example was given whereby I^Tj^Xcode 
was automatically generated to draw the diagrams. A suggested further application 
is automatic code generation to evaluate the series in specific cases. 

Acknowledgements 



*available on request from the author 



Automatic Generation of Vacuum Amplitude 



This work was supported by the UK EPSRC 

Appendix A 

module hug 
type pair 

integer : : vl , v2 
end type pair 
type diagram 

integer, dimension( : ) , pointer :: occ 
end type diagram 
integer :: npairs, norder 

contains 

logical function consistent (diag,nlines .pairs) 
implicit none 

type (diagram) , intent (in) :: diag 

integer, intent (in) :: nlines 

type(pair), dimension( : ) , intent(in) :: pairs 

integer :: linesum, i, vertex, lines, j, start, fin 

integer, allocatable, dimension(:) :: link 

logical, allocatable, dimension(:) :: linked, to_try, tried 
logical : : isdiaglinked 

consistent = .true. 

! check each vertex for correct number of lines: 

allocate (link(norder) ) 

link=0 

do i=l , size(diag7 occ ( : ) ) 

link (pairs (i)7 vl)=link (pairs (i)°/ovl)+diag / occ(i) 
link (pairs (i)7.v2)=link (pairs (i)°/,v2)+diag°/,occ(i) 

end do 

do i=l, norder 

if (link(i)/=4) then 

consistent=. false. ; return 

end if 
end do 

deallocate (link) 

! Test for unlinked diagrams. 

isdiaglinked= . false . 

allocate (linked (norder) ,to_try (norder) , tried (norder) ) 
linked=. false. ; to_try= . false . ; tried= . false . 
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linked(l)=.true. ; to_try(l)=.true. 
do 

vertex=0 

do i=norder,l,-l 

if (to_try(i) .eqv. .true.) vertex=i 
end do 

if (vertex==0) exit 
do i=l , size (pairs) 

if (diag7,occ(i)/=0) then 

if (pairs (i)7«vl==vertex) then 
linked(pairs(i)7 v2) = .true, 
if (tried (pairs (i)7 v2) .eqv. .false.) then 

to_try (pairs (i)°/ v2) = .true, 
end if 

else if (pairs (i)°/ v2==vertex) then 
linked(pairs(i)°/ vl) = .true, 
if (tried (pairs (i)°/ vl) .eqv. .false.) then 

to_try (pairs (i)°/.vl) = .true, 
end if 
end if 
end if 
end do 

tried(vertex)= .true . 
to_try (vertex) = .false . 
isdiaglinked= . true . 
do i=l,norder 

if (linked(i) .eqv. .false.) isdiaglinked= . false . 
end do 
end do 

deallocate (linked, to_try, tried) 
if (isdiaglinked. eqv. .false.) consistent= . false . 
end function consistent 

subroutine label (diag, pairs) 
implicit none 
type (diagram) : : diag 

type(pair), dimension( : ) , intent(in) :: pairs 
integer, dimension( : ) , allocatable :: nups 
integer : : offset = 

allocate (nups (npairs) ) 

call newlabel (of f set ,nups ,diag°/ occ ( : ) , pairs) 
deallocate (nups) 
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end subroutine label 

logical function testlab(nuparray , pairs, occarray) 
implicit none 

type(pair), dimension( : ) , intent(in) :: pairs 

integer, dimension( : ) , intent(in) :: nuparray, occarray 

integer :: vstart, vend, i 

integer, dimension( : ) , allocatable :: enter 

allocate (enter (norder) ) 

enter=0 

do i=l,npairs 

vstart = pairs(i)°/„vl 
vend = pairs (i)7 v2 
if (vstart < vend) then 

enter(vend) = enter (vend)+nuparray(i) 
enter (vstart) = enter (vstart) + occarray (i) -nuparray (i) 
else 

enter (vstart) = enter (vstart)+nuparray(i) 
enter(vend) = enter(vend)+ occarray(i)-nuparray(i) 
end if 
end do 

testlab= .false . 

if (all (enter==2) ) testlab = .true, 
deallocate (enter) 
end function testlab 

recursive subroutine newlabel (of f set , nuparray, occarray, pairs) 
implicit none 
integer : : offset 

integer, dimension(:) :: nuparray, occarray 

integer :: nuplo, nuphi , n, j 

type(pair), dimension( : ) , intent(in) :: pairs 

if (of f set==npairs) then 

if (testlab (nuparray , pairs, occarray) . eqv . .true . ) then 

write (*, ' (al,lx,100(il)) ') '*', (nuparray (j ), j=l ,npairs) 

end if 

return 
end if 

of f set=of f set+1 

nuplo = max (occarray (off set) -2,0) 
nuphi = min(occarray (of f set) ,2) 
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do n = miplo,nuphi 

nuparray (of f set) = n 

call newlabel (of f set, nuparray, occarray, pairs) 
end do 

of f set=of f set-1 
end subroutine newlabel 

recursive subroutine new(nlines, offset, temp, pairs) 
implicit none 

integer, intent(in) :: nlines 

type (diagram) : : temp 

integer, intent (inout) :: offset 

type(pair), dimension( : ) , intent(in) :: pairs 

integer :: i, maxlinks =3, j, ilo = 

if (of f set==nlines/2-l . and. sum(temp7 occ (1 : offset) ) . ne .4) return 
if (of f set==npairs) then 

! We have got to the end of the array - check for consitency: 
if (consistent (temp , nlines , pairs ) ) then 

write (*, ' (al,lx,100(il)) ') ' + ' , (temp7.occ(j) , j = l,npairs) 
call label (temp, pairs) 
end if 
return 
end if 

if (npairs==l) maxlinks = 4 
do i=ilo , maxlinks 

temp7 occ (1+of f set)=i 

if (sum(temp°/„occ (1 : of f set+1) ) >nlines) cycle 
of f set=of f set+1 

call new(nlines, offset, temp, pairs) 

of fset=off set-1 
end do 
end subroutine new 
end module hug 

program hugenholtz 
use hug 
implicit none 

integer :: nfixed, nlines, index, i, j, offset = 
type(pair), dimension( : ) , allocatable :: pairs 
type (diagram) : : temp 

write (unit=*,fmt= , ("Input order of diagrams: ") ' , advance= 'no ' ) 
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read(unit=* ,fmt=*) norder 

npairs = norder * (norder - 1) / 2 

allocate (pairs ( 1 : npairs) , temp7 occ ( 1 : npairs) ) 

nlines = 2 * norder 

index = 1 

do i = 1 , norder 

do j = i + 1 , norder 

pairs (index) 7 vl = i ; pairs ( index) 7 v2 = j 
index = index + 1 

end do 
end do 

tempZocc ( : ) = 

call new(nlines, offset, temp, pairs) 
deallocate (temp %occ .pairs) 
end program rmgenholtz 
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Figure 3: The 39 fourth order diagrams 



